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Abstract. - We present numerical and analytical results describing the effect of hydrodynamic interactions on the 
dynamics of a short polymer chain in solution. A molecular dynamics algorithm for the polymer is coupled to a direct 
simulation Monte Carlo algorithm for the solvent. We give an explicit expression for the velocity autocorrelation 
function of the centre of mass of the polymer which agrees well with numerical results if Brownian dynamics, 
hydrodynamic correlations and sound wave scattering are included. 



Introduction. - The aim of this paper is to present numerical and analytical results describing the effect of 
hydrodynamic interactions on the dynamics of a short polymer chain in solution. Modelling a dilute polymer 
solution is a difficult task. This is because the dynamical behaviour of the polymer can be dominated by hydrody- 
namic interactions between different parts of the polymer chain. These interactions are long-ranged and develop 
very slowly compared to the time scale of the Brownian fluctuations of individual monomers. Thus if a molecular 
dynamics time-step is chosen to correctly follow the Brownian dynamics it becomes prohibitively expensive in 
computer time to also model the dynamics of enough solvent molecules for sufficient time to allow hydrodynamic 
correlations to develop. 

To move towards overcoming this problem we introduce a hybrid numerical approach. The dynamics of the 
chain is treated exactly (within the limitations of a finite time step) by a molecular dynamics solution of Newton's 
equations of motion. The solvent is modelled using a direct simulation Monte Carlo algorithm [Q . The algorithm 
is constructed in such a way that energy and momentum are conserved and that at equilibrium the system is 
described by a microcanonical distribution. However, molecular details of the solvent are not included: it acts as 
a momentum-conserving heat bath which can support hydrodynamic modes. 

The physical quantity we choose to measure in the simulations is the velocity autocorrelation function of the 
centre of mass motion of a polymer chain. This provides clear evidence of a fast exponential decay at early times 
which results from collisions in which the solvent particles are uncorrelated and a slow algebraic decay at later 
times arising from hydrodynamic interactions that develop in the fluid. We propose a theory which reproduces 
the behaviour of the velocity autocorrelation function and demonstrate that both limits arise naturally from a 
simple equation of motion which incorporates the solvent-polymer interaction as a viscous drag term. We also 
argue that, as the fluid is compressible, a feature common to many mesoscale approaches, both transverse and 
longitudinal velocity modes participate in momentum transfer yielding a significant correction to the value of the 
chain diffusion coefficient p|. The theoretical prediction is shown to agree well with the numerical results. 

The paper starts with a description of the numerical method. Next we present an analytic approach to the 
calculation of the velocity autocorrelation function. We emphasise that both the Brownian and the hydrodynamic 
contributions are important, and that the latter includes contributions from sound waves. We show how the 
Kirkwood formula for the diffusion constant of the centre of mass of the polymer chain follows from the theoretical 
development and hence discuss the relation of our work to previous mesoscale simulations of polymer dynamics 
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Direct simulation Monte- Carlo algorithm. - We construct a hybrid algorithm for the dynamics of polymers 
in solution by exploiting the time-scale separation between the microscopic motion of the polymer beads and the 
much slower propagation of hydrodynamic modes within the fluid. The evolution of the polymer beads follows 
a molecular dynamics algorithm. The solvent is modelled on a mesoscopic length scale using a direct simulation 
Monte Carlo approach. 

Consider first the solvent. The direct simulation Monte Carlo algorithm is based upon the motion of N s 
particles in continuous space but discrete time. Let the particles, labelled i = 1, 2 . . . N s , take positions Xi(t) with 
velocity Ui(t) at time t. The evolution for unit time comprises a free streaming step applied to each particle i 

X 4 (t+l)=3C i (t)+U i (t) (1) 

followed by a collision step. The collision step must (i) preserve mass, momentum and energy locally so that the 
macroscopic equations of motion are obtained in the continuum limit; (ii) allow momentum transfer between the 
particles so that equilibrium can be achieved. Many different schemes are possible. Here we follow Malevanets 
and Kapral j^] and divide space into L x x L y x L z cubic cells each containing, on average, N s / L x L y L z particles. 
The secular velocities of particles in each cell are then rotated according to 

Ui(f + 1) = U(i) + ft(ui(i) - U(t)) , (2) 

where U is the velocity of centre of mass of particles in each cell and 1Z is a rotation matrix. In the present study 
1Z is taken to be a rotation by 7r/2 around an arbitrarily chosen axis. 

Various hydrodynamic properties of the solvent are needed to compare the analytical and numerical results. 
These depend on the details of the collision step and are most easily obtained numerically using either the Green- 
Kubo formalism or non-equilibrium dynamics simulations. Parameters for the solvent component are temperature 
kgT = 1.0, mass m = 4.0 and density p = 2.0. A time integral over the stress autocorrelation function gives the 
kinematic viscosity coefficient as v = 0.25. The decay rate of a small longitudinal sinusoidal perturbation gives 
the attenuation rate of a sound wave to be T = 0.275 and the speed of sound c = 0.645. 

A chain is introduced into the system by choosing N beads which are connected by the harmonic potential 
W n — K,\\r n — r n+ i|| 2 /2. Apart from the harmonic potential the polymer beads are taken to be non-interacting to 
facilitate comparison to the predictions of the Zimm and Rouse theories . The beads are included in the same 
collision step (^) as the solvent particles. In between the collision steps the polymer dynamics is integrated using 
a velocity Verlet algorithm H with time step At. 




t 



Fig. 1 - Velocity autocorrelation function of a polymer chain. The circles show the computed values and the solid line 
represents the theoretical curve calculated using equation ([]). 



The simulations were used to measure the velocity autocorrelation function of the centre of mass of the polymer 
chain at equilibrium. The results are presented in Fig. [I]. The parameters used in simulations were AI = 8.0, 

N = 128 and k = 3.0. The radius of gyration of the chain was R g = J = 2.31. 
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Calculation of the centre of mass velocity autocorrelation function. - The aim of this section is to derive a 
formula for the velocity autocorrelation function of the centre of mass of the polymer chain. We commence by 
writing the equation of motion of the monomers as [bl 



M— v„(t) = f n (t) + f(u(r n (t), t) - v n (f)) 



(3) 



where v„ is the velocity of nth monomer and r„ its position at time t. f„ is the force on the nth monomer exerted 
by the other beads in the chain. u(r„(t),£) is the velocity of the solvent at r n at time t. The coefficient £ gauges 
the strength of viscous friction between the chain and the solvent. 

We are interested in the behaviour of the velocity of the centre of mass of the polymer chain V(t). Summing 
(||) over all beads, we note that the contribution from the internal forces cancel. The resulting equation can be 
integrated to give 



V(i)= e ^/ A 'V(0) + JL J dTe -t(t-r)/MJ2u(r n (T),T) 



(4) 



The next task is to calculate u(t w (t), r). This can be done by using the linearized thermo-hydrodynamic 
equations for a compressible fluid [Q. We assume the transfer of momentum from the chain to the fluid is 
sufficiently fast that 



u(x,0) = — y>„<5(x-r n (0)) 
P 



(5) 



Evolution of the fluid velocity field is then given by [[HI 



. ^ M 

la[K,t) = } 

„ P 



-j^(k 2 5 al3 - k a k )e vkH + ^k a k p e I7 ''' cos /,-<•/ 



»ik.r„(0) 



«n/s(0) 



(G) 



Using the identity 



u(r„(i),t) 



(27T) 



^3 J dke- lk r "Wu(k,t) , 



we can now calculate the chain centre of mass velocity autocorrelation function 

C{t) = I (V(t)V(0)) , 



(7) 



(8) 



where the angular brackets denote an equilibrium ensemble average over the Gibbs distribution. Substituting the 
expression for the fluid velocity (Q) into the formula for the polymer centre of mass velocity (||) and using the 
result in rtq) gives 



C(t) 



-£t/M 



MN P M{2tt) 3 N 



dre - m -r)/M / dk 



-vk A T 



~ rfc2r cos kcT 



5(k,r) 



(9) 



where we have used the usual definition of the structure factor 



,ik.(r n (T)-n(0)) 



(10) 



We assume that the structure factor exhibits a two-stage decay; an initial fast decay due to the diffusion of 
individual beads followed by a later stage when it decays as exp(— Dk 2 t). We found that in our simulations the 
structure factor is well approximated by 



S(k,t) = Ne -Dk 2 t e ~R g k*/3 



(11) 
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where D <C v. This expression is obtained assuming that the chain configurations at times and t are uncorrelated. 

The behaviour of C(t) obtained from a numerical integration of equation (^) using (|ll]) is compared to the 
simulation results in Fig. |l|. Agreement is good. Small deviations of the theoretical curve from the numerical 
results on short-times scales can most probably be attributed to a breakdown of the assumption of a single friction 
coefficient with no memory effects in equation (^). 

It is sometimes helpful to divide C(t) into a short time or Brownian and a long time or hydrodynamic 
contribution. These can be identified as the first and second terms in equation (^). Note that there is no 
clear time scale separation between the two contributions: we expect this to be the case for the short polymers 
considered in mesoscale simulations. Note also that both contributions follow from the same physical mechanism, 
the viscous drag between the fluid and the polymer. However, the first term corresponds to Brownian motion when 
the fluid molecules can be considered to move at random; the second occurs when hydrodynamic correlations in 
the fluid lead to interactions between the monomers. It can be "switched off" in the simulations by randomizing 
the solvent molecules' velocities after each collision step. 

The second term in the square brackets in (Q) arises from sound wave-chain interactions. This coupling may 
lead to a negative minimum in the velocity autocorrelation function Q . 

Calculations of the chain centre of mass diffusion coefficient. - In this section we estimate contributions from 
different terms of equation (^) to the diffusion coefficient of the centre of mass of the polymer chain. Replacing 
S(k, t) in equation IM) by its value at some intermediate time <5'(k, t 1 ) and integrating gives 



D 



k B T 



1 



67:77 \ \ ri (t) -r n (t')\ 



k B T / e - C |r,W-r„(t')|/r' 
|r,(t)-r„(i')| 



12irpT \ 



(12) 



In equation (|l2j) the first term corresponds to the microscopic or uncorrelated contribution and has the same 
scaling as Rouse's model predicts. The second term gives the Zimm contribution due to hydrodynamic interactions 
between beads and the last term arises from scattering of sound waves on the chain as the longitudinal component 
of momentum drains from the system. 

For the simulations the first and third terms in this equation are approximately 10% of the total value. The 
same order corrections are expected in water, where typical parameters are c ~ 10 3 ms _1 and T ~ 10~ 6 m s , 
for chains with radius of gyration R g ~ 10~ 8 m. For longer chains when the condition R g c ^ T is satisfied the 
sound contribution to the diffusion coefficient vanishes. 
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Fig. 2 - Diffusion coefficient as a function of system size. Solid circles correspond to numerical values of the diffusion 
coefficient for systems of size 8, 16 and 32. The line is a fit to the data. 



Hydrodynamic interactions are long range and finite-size effects lead to large corrections to the diffusion 
coefficient ]Tl| ] . In order to estimate this effect we replace integration over k in equation (Q) with the summation 
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over discrete wave-vectors k = (27r/L)n, k ^ 0. Then in the intermediate time regime M/£ « i « R 2 g l v we 
obtain 

c *® = c ®-jI$Fn I dk l siKt) (13) 

\k\i<ir/L 

Hence, assuming that the structure factor varies only weakly with time, 

1 2.44 



2 k B T f 1 k B T 



3 (27r) 3 77 J k 2 Q-kt] 

\k\i<n/L 



Rh 



(14) 



This form of system size dependence is consistent with the correction obtained by Dunweg et al. jll[] using different 
arguments. In Fig. ^ we present values of the diffusion coefficient of a chain obtained in simulations of systems 
with dimensions 8 3 , 16 3 and 32 3 . A numerical fit gives for the finite size correction 

where D = 1.05 x 10~ 2 . Alternately we may compute the hydrodynamic contribution to the diffusion coefficient 
from equation (|T^ ) . Averaging yields 

kfiT o 
D= J= — = 1.15 x 1CT 2 . (16) 



This is consistent with the numerical value of the diffusion coefficient obtained in the simulations within statistical 
errors and the difficulties inherent in calculating the finite-size corrections. 

Discussion. - We have described a new mesoscale method to simulate the dynamics of polymer chains in 
solution. The polymer chain itself is modelled using standard molecular dynamics. This has the advantage that 
the architecture of the chain and the interactions within it can be varied. The solvent is described by a direct 
simulation Monte Carlo algorithm which reproduces the thermohydrodynamic equations of motion on sufficiently 
long length scales. By ignoring the molecular detail of the solvent it becomes feasible to study the build up of 
hydrodynamic correlations between different monomers. 

The numerical results were shown to agree well with an equation of motion for the monomers which assumed a 
linear coupling to the velocity of the surrounding fluid via a friction coefficient £. The velocity-velocity correlation 
function comprises two contributions; an initial exponential decay resulting from Brownian collisions of monomers 
with uncorrelated solvent molecules and, at longer times, a hydrodynamic contribution as monomers interact via 
hydrodynamic modes propogated through the solvent. There is no clear time-scale separation between the two 
contributions for polymer chain lengths accessible to computer simulations. Hence the Zimm limit, which assumes 
that the velocity-velocity correlation function has decayed to zero is not appropriate for the systems we consider 
here. However a Kirkwood formula for the diffusion constant (jl^) follows naturally and shows that the Zimm 
scaling D ~ Rq 1 for the hydrodynamic contribution to the diffusion constant is still expected. 

Other mesoscale approaches which have been used to study dilute polymer solutions include dissipative particle 
dynamics |^,[5| and a polymer chain updated using molecular dynamics coupled to a lattice Boltzmann solvent 
Both approaches gave good results for the Zimm scaling of the diffusion constant. Qualitative values were 
consistent but precision was limited by the difficulty in separating out the Brownian contribution from the 
hydrodynamic one and because of strong finite-size corrections. It would be of interest to measure velocity- 
velocity correlation functions using these methods to enable a more detailed comparision of the dynamics to that 
obtained using the mesoscale approach described here. 

Several open questions remain. We have concentrated on the dynamics of the centre of mass of the polymer 
chain. It would be interesting to next consider the internal dynamics (Rouse modes) which will be affected by 
the interaction between the polymer and the solvent on shorter length and time scales. It will also be possible to 
investigate polymer dynamics under shear and the effect of the architecture of the polymer chains on their flow 
behaviour. 
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